Controlling One-Dimensional Langevin Dynamics on the Lattice 
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Stochastic evolutions of classical field theories have recently become popular in the study of prob- 
lems such as determination of the rates of topological transitions and the statistical mechanics of 
nonlinear coherent structures. To obtain high precision results from numerical calculations, a careful 
accounting of spacetime discreteness effects is essential, as well as the development of schemes to 
systematically improve convergence to the continuum. With a kink-bearing (j> 4 field theory as the 
application arena, we present such an analysis for a 1+1-dimensional Langevin system. Analyt- 
ical predictions and results from high resolution numerical solutions are found to be in excellent 
agreement. 
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I. INTRODUCTION 

In recent years there has been growing interest in ex- 
tracting non-perturbative quantum dynamical informa- 
tion such as topological transition rates from numerical 
Langevin and Monte Carlo solutions of classical field the- 
ories at finite temperature [Q. At the next level of so- 
phistication, several attempts have been made at devel- 
oping schemes that treat low-lying modes classically and 
high frequency modes quantum mechanically ||. More- 
over, the equilibrium and nonequilibrium classical sta- 
tistical mechanics of nonlinear coherent structures such 
as kinks has historically received much attention || in 
the condensed matter literature. Until fairly recently, 
computer memory and performance restrictions were suf- 
ficiently severe that Langevin evolutions could only be 
carried out at fairly low levels of accuracy and resolu- 
tion. However, present-day supercomputers have over- 
come this problem, at least for low dimensional problems, 
and one can well contemplate systematically studying, 
understanding, and improving the accuracy of stochastic 
evolutions. In this paper we present just such a study 
applied to 1+1-dimensional Langevin systems. 

Our focus will be on lattice errors for quantities com- 
puted at thermal equilibrium. In calculations of this 
type, a stochastic partial differential equation (SPDE) 
augmented with a fluctuation-dissipation relation is 
solved as an initial value problem using finite differences. 
Because of the fluctuation-dissipation relation, the sys- 
tem is eventually driven to thermal equilibrium and at 
late times one may measure values of thermodynamic 
quantities as well as time and space dependent correla- 
tion functions. These quantities can depend on the lattice 
spacing, Ax, on the total system size, on the discretiza- 
tion used for spatial operators, and on the timestep- 
ping algorithm used to solve the resulting set of coupled 
stochastic ordinary differential equations. In one space 
dimension a fairly complete description can be given since 
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the question of lattice effects is one of convergence prop- 
erties of SPDEs rather than of renormalization. 

The configurational part of the partition function of a 
classical field theory in one space dimension is free from 
divergences. In particular, quantities such as kink densi- 
ties, measured from finite difference solutions of the cor- 
responding SPDEs, converge to a well-defined limit as 
the lattice spacing is reduced towards zero. The question 
of exactly how the convergence scales with Ax is still a 
matter of practical importance: numerical solutions are 
limited by the available computing power and memory to 
a finite range of values of Ax. While finite volume effects 
can be important in small lattices they are not important 
if the lattice size is much larger than the longest correla- 
tion length. We will assume that this is always the case 
in the considerations below. 

A complete constructive procedure for determining the 
spatial lattice error, and possibly eliminating it to some 
order in Ax, exists. The method proceeds as follows. 
In equilibrium, the probability of a given set of config- 
urations can be calculated from the static solution of 
the Fokker-Planck equation corresponding to the particu- 
lar spatial discretization and time-stepping algorithm ap- 
plied to the SPDE of interest. With time-stepping errors 
tuned to be sub-dominant, the transfer integral Q cor- 
responding to the lattice Hamiltonian can be evaluated 
to some given order in Ax. Correlation functions and 
thermodynamic quantities, which can all be extracted 
from the transfer integral, explicitly exhibit lattice de- 
pendences allowing the convergence to the continuum to 
be read off directly. We describe this procedure in more 
detail below. 

Discreteness effects have been considered before in the 
context of kink dynamics ||. Trullinger and Sasaki 
have already obtained the lowest-order discreteness cor- 
rections to the Schrodinger equation that emerges from 
the transfer integral approach |^. They found that the 
lowest order correction is of order Ax 2 and, in first-order 
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perturbation theory, is equivalent to a corrected effec- 
tive potential. As we show below, the latter result can 
be adapted not only to compute the order of the lat- 
tice errors but also to introduce a local counterterm in 
the stochastic evolution equations that can drastically 
improve the convergence to the continuum. Our results 
from high resolution numerical solutions are in excellent 
agreement with the theoretical predictions. 

The class of problems considered here are 1+1- 
dimensional classical field theories defined by the Hamil- 
tonian: 



H 



dx 



K 2 



1 /<9$ 



2 U) 



(1) 



The corresponding continuum SPDE 



(2) 

is second order in time, where with = 1/kT, the noise 
and damping obey a fluctuation-dissipation relation: 



(F(x, t)F(y, a)) = 2 V /3- 1 S(t - s)S(x - y) 



(3) 



In this paper we will adopt the example of the double- 
well $ 4 theory: V($) = -(m 2 /2)$ 2 + ( 5 2 /4)$ 4 . We shall 
work in a dimcnsionlcss form of the theory given by the 
transformations: = &/a, x = mx, and t = mt, where 
a 2 = m 2 /g 2 . Under these transformations, the original 
Hamiltonian becomes H = H/(ma 2 ) where H is of the 
same form as the original Hamiltonian H, except that 
the potential V{<j>) = -(l/2)0 2 + (l/4)0 4 . 

This theory admits the well-known (anti-)kink solu- 
tions which, at zero temperature, are exact solutions 
of the static field equations connecting <j> = — 1 at 
x = (+) — co to <p = +1 at x = {— ) + oo. In ther- 
mal equilibrium, the balance between noise and damping 
is manifested in the balance of nucleation and annihila- 
tion of kink-antikink pairs @). At low temperature, WKB 
techniques applied to the transfer integral |||| yield the 
following approximation for the density of kinks: 



Pk ^ (Ek/kT) 1 / 2 eM-Ek/kT), 



(4) 



where Ek = \/8/9, the energy of an isolated kink. Sup- 
porting numerical evidence exists [ flO| , but precise results 
have been difficult to obtain until recently due to the 
large amount of computing time needed at temperatures 
low enough to clearly distinguish kinks. The best results 
obtained so far are for a special double-well potential 
where exact theoretical computations can also be carried 
out. In this case, it has been shown that the theoretical 
and numerical results agree within statistical bounds set 
by the finite volume of the simulations ITT). 

The classical partition function for a <jr theory, in any 
spatial dimension 2 < D < 4, is super-renormalizablc, 



i.e., there are a finite number of perturbative diagrams 
that are divergent in the continuum, but can be appro- 
priately subtracted by the inclusion in the theory of a 
finite number of suitable counterterms. The situation is 
different for D = 1: the continuum partition function is 
finite and no renormalization is necessary. 

An alternative approach to the one described here has 
been suggested by Gleiser and Miiller jl2| who have 
proposed a perturbative counterterm for use in 1+1- 
dimensional Langevin equations. A weakness of the lat- 
ter proposal is that it relics on an approximation to the 
free energy; in many situations the latter is a poor indica- 
tor of the true dynamics of field theories jlj|. Moreover, 
their counterterm is based on an approximate effective 
potential calculated by perturbing about a uniform state. 
We will show below, with both analytic and numerical 
results, the inadequacy of perturbative counterterms in 
dealing with the convergence to the continuum. 

The paper is organized as follows. In Section II we 
consider the evolution of the probability density of the 
discretized SPDE. We use the method of Horowitz |TJJ] 
to examine the effect of time discretization on the equi- 
librium density. The transfer integral is introduced in 
Section III. We perform calculations at finite Ax and 
show that the leading order corrections to the contin- 
uum of observable quantities are proportional to Ax 2 . 
Examination of the form of the Schrodinger equation at 
finite Ax reveals a natural choice for a local countert- 
erm with which to improve the convergence properties 
of discretized Langevin equations. The alternative one 
loop approach of Gleiser and Miiller is examined in Sec- 
tion IV. Numerical results are presented in Section V. In 
Section VI we end with a discussion of our results. 



II. THE DISCRETE TIME FOKKER-PLANCK 
EQUATION 

Our first step in determining the (equilibrium canon- 
ical) distribution to which a given Langevin dynamics 
converges for long times is to derive the corresponding 
Fokker-Planck equation. This can be done on the lattice 
as well as in the continuum. 

On the lattice, an SPDE is solved by updating 2N 
quantities {4>i{t), ni(t)} where i = 1, . . . , N. We take the 
lattice Hamiltonian in one space dimension, £/i a t, to be 



N 



with 



Hi 



lal 



= AxJ2 



\A + s(<t>i) 



1 (<t>i+l ~ 4>i) 

2 Ax 2 

-^<Pi + -7<Pi ■ 



+ V(<f>i), 



(5) 

(6) 
(7) 



The corresponding Fokker-Planck equation for the 2N 
variables has a static solution that can in principle be 
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attained at late times in a Langevin simulation (in the 
sense of ensemble averages over individual simulations) . 

In practice the time as well as the space discretization 
of a Langevin equation leads to errors. The simplest 
stochastic timestepping is of the Euler type and can be 
written as: 



%i(t + At)= m(t) - At 



dH ut 



dMt) 



6(f), 



<f>i(t + At)= + AtTT.it). 



(8) 



We have chosen the case of additive Gaussian white noise, 
related to the damping 77 by the (suitably discretized) 
fluctuation-dissipation relation: 



&(*)> = 0, <&(t)&(f)> = 



2t] 1 
~P At Ax 



-Sin5 



ijOtt'- 



(9) 



In order to understand the effect of time discretization, 
it is possible write a discrete time Fokker-Planck equa- 
tion, describing the evolution of the probability density 
functional associated with (||)-(^|) Jl^ |: 



P[{k, 4>},t- 
d 



At] = exp -At 



d dH, 



lat 



84>i diTi 



exp 



At 



lat 



lat 



. r] d 2 
+ Atj-—^ 



(10) 

P[{ir,4>M, 



where summation over repeated indices is implied. For 
simplicity this will be assumed in what follows and in- 
dices dropped. The discrete time equation ( |l0| ) can be 
written in the form 



F[{vr,0},t + At] = e 



_ -AtH FP 



P[{n, </>},*]• 



(11) 



The operators in the two exponents in (|i0| ) are non- 
commuting. To reduce (|l^) to the form (^L 1|) we use the 
Campbell-Bakcr-Hausdorff theorem: given the operators 
A and B, there is (formally) an operator C such that 



e A e B = e u , with 

C= A + B + ±[A 1 B] + ±[A,{A,B}} 
+ L[B,[B,A]] + ... 
Expanding to first order in At, we have 



(12) 
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dir da 



0_ 

On 



lat 



lat 



dHiat r) d 
71 ^ + Pd^ 



dir 
dH ut 



d(p 



dTT 

dH lat d 2 H 



(13) 



lat 



d 



dir dej) 2 dir 



0(At 2 



Notice that each factor of H\ at introduces a power of Ax. 

The solution of HppP[{(f>, ir}} = is the canonical dis- 
tribution approached by the discretized system at late 



times. Its form can be computed for small At. To ze- 
roth order for the momenta and order At for the fields 
we obtain 

F[{tt, 0}] = exp ( - Ax ]T + (3S(&) (14) 



B dS 



where /?' = 0(1 + Ati). Note that the discretization 
induces cross terms between <fi and n in the canonical 
distribution. This is a general feature of higher order 
solutions in At. (These terms rapidly become very com- 
plicated.) Different time discretizations lead to different 
discrete time Fokker-Planck equations. The numerical 
simulations described below employed a stochastic sec- 
ond order Runge-Kutta algorithm 

The equilibrium density of configurations of the space- 
and time-discretized theory, is obtained by performing 
the Gaussian integral over the momenta in (|14|): 



P[{<f>}]= A exp 



-f3AxJ2 




(15) 



The effect of the time discretization is explicitly seen as a 
modification of the equilibrium density. Further integra- 
tion of ( |l4| ) cannot be performed so easily because each 
S(<fii) depends also on 4>i+\. 

The functional integral of P[<f>] over (f> defines the con- 
figurational partition function 



/ D4>e 



-ps[<f>] 



N 

n 

i=l 



(16) 



which we study in the next section. Here &<f>i — Ndfa, 
with N = 



2ttAx ' 

In principle, the partition function as calculated above 
would include artifacts from both the time and space 
discretizations. In actual computational practice, given 
that we can estimate the order of the time stepping errors 
as shown already, it is not difficult to reduce the time 
step to a level where the remaining errors are suppressed 
compared to the errors from the spatial discretization. 
Once this is done, we may safely ignore the discretization 
in time and concentrate solely on the errors due to the 
spatial lattice. 



III. THE TRANSFER INTEGRAL 

To compute the partition function explicitly we make 
use of the transfer integral method Q. The configura- 
tional partition function is given by 



Z,i 



/ d0i...d0jvIjT(0i,0,' + i), 

J — OO • 1 



(17) 
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where 

T((/>i,<j> i+ i) = 
[ 1 



exp 



-8 Ax 



H+i 



Aa 



+ V(4h) + V(<Ih+i) 



and 4*n+i — 4>i implements spatially periodic boundary 
conditions. The difficulty with evaluating Z^, lies in the 
coupling of integrals at different space points. The idea 
behind the transfer operator method is to "localize" the 
evaluation of the integrals in (|l7|). 

The transfer operator T is defined as follows 



fip(4> i+1 ) = / d^iT(<f>i,(f)i + i)tp(<f)i). 



(18) 



Suppose we can find the eigenvalues of T. That is, sup- 
pose we can solve the following Fredholm equation: 



d4>iT(^i,(f)i + i)i) n (4>i) = t n il} n (4> i+ i) , (19) 



where the t n are positive constants that we write for later 
convenience as 



Then 



EC 



(20) 



(21) 



In the limit N — > oo, the sum (^lj) is dominated by the 
largest eigenvalue: 



/3Le 



(22) 



where L = NAx is the physical length of the lattice. In 
the thermodynamic limit L — + oo, the free energy density 
is simply = £q. It is clear that once the partition func- 
tion is known in the thermodynamic limit we may com- 
pute from it any thermodynamic quantity. Moreover, it 
is possible to show that spatial, and in linear response 
theory, temporal correlation functions can also be com- 
puted via a knowledge of the spectrum of the transfer 
operator flili|| . 

We now turn to the procedure for solution of ( |l9| ) by 
first converting it into an infinite order partial differential 
equation. We first rewrite (Q9h as 



e 2 



e -^w. + i-wy^-0> + i)dA?0 1+lx( ^ +i) 



-[3Ax£ n 



VVi(0 



where 



X{<t>) 



(23) 

(24) 



The special form of the Fredholm kernel has led to a 
simple Gaussian integral that yields: 



-i/3AxV(0 i+1 ) (Ax/20)d 2 /d</> 2 i+1 



(e s 



0AxV(<p i+1 ) 



— e -PAxe r 



ll>n(<t>i+l) 



(25) 



This (exact) result yields the form e e e ip — e ip 
where U and D are operators and C is a real number. 
The Campbell-Bakcr-Hausdorff scries in this case is for- 
mally an expansion in powers of Ax. To linear order in 
Ax, the CBH expansion applied to (123) yields: 



-f3AxV(<t>)+Ag.a 

e ^^WnW) 



-PAxtn 



or equivalently 
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(26) 



(27) 



The transfer integral technique thus reduces the calcula- 
tion of to the calculation of the eigenvalues e„ of a 
corresponding Schrodinger equation: 



where {7(0, 0) = V(<p). The calculation is explicitly per- 
formed on the lattice, at finite Ax: leading order correc- 
tions to the eigenvalues of the Schrodinger equation ( g^ ) 
are proportional to Ax 2 . For the problem at hand, one 
finds [HI, 



+ V(<P) + ^(Ax) 



dV 



i d 2 

28 2 dcj) 2 

1 d 2 V d 2 
'2{P 7302" • 28 2 d(f> 3 d(f> ' 88 2 d(f> 4 



1 d 3 V d 1 d A V 
+ 



(29) 



ip n = e n ip n 



Higher order corrections in Ax in (|29|) can be computed 
in a tedious though straightforward fashion by going to 
higher orders in the CBH expansion. It is easy to show 
from the symmetric form of (|2^) and the Hermitian/anti- 
Hermitian alternation of terms in the CBH expansion 
that the error terms are always of even order in powers 
of Ax. Thus, if a method is found to cancel errors up to 
a certain order m, it automatically reduces the error to 
order m + 2. 

The simplest example of Ax dependence is the free 
theory: V = \cj) 2 . Then @ reduces to 



1 d 2 

2/3' 2 d(j> 2 



1 

—m cp 



(30) 



with 8' = 8/VT^, m! = y/\ + a/2 and a = (Ax) 2 /6. 
This implies in particular for the energy spectrum e n — 

( n + 1) f- - ( n + i 1 - M Ax ) 2 )- The free theor y 

is a convenient special case because the corresponding 
SPDE is linear and exactly solvable. Quantities such as 



4 



(</> 2 (cc)) can be evaluated exactly and compared to the 
results from the transfer integral. Both procedures agree, 
and 



{<?{x)) = 



1 



1 



2ft 



Ax 2 
4 



1 . , 

-Ax' 



0(Ax 4 



(31) 



Note that the leading dependence on lattice spacing is 
proportional to Ax 2 . 

We now turn to the question of lattice errors in de- 
termining the kink density, which, at sufficiently low 
temperatures, is controlled completely by the correlation 
length derived from the two-point function (0(0) 4>{x)). 
Applying the transfer integral formalism, it is easy to 
show that this correlation function is a sum of exponen- 
tials with exponents proportional to differences of eigen- 
values of (p8|). The correlation length is determined by 
the energy difference between the ground and first ex- 
cited states of (p8| ) ||. At low temperatures, the WKB 
(or semiclassical) approximation is excellent and this en- 
ergy difference is the exponentially small tunnel-splitting 
term. Note that at low temperatures the kink density is 
|iven directly by the correlation length, pk — l/(4Aoo) 



At low temperatures the first two eigenf unctions of (|27 
are of the form 



(32) 



where ips is the (symmetric) ground state and tpA is the 
(antisymmetric) first excited state. Here ipL and ipR are 
the usual localized states, one on each side of the barrier. 
To estimate the error due to finite lattice spacing we use 
first order perturbation theory in (Ace) 2 . The corrected 
energies are then, 



E 
E 



(Ax) _ 



(Ax) 



= Eo + (1> 8 \8H\ilis), 
Ex + (1> a \6E\iI>a), 



(33) 



where E$ and E\ are the results from the continuum 
theory and SH ~ (9(Ax 2 ) is the error Hamiltonian. It 
follows that the energy difference is 



AE 



(Ax) 
10 



AE 10 -2(<iP l \SH\tPr). 



(34) 



The error Hamiltonian can be read off from (29) and it is 



clear that the error in energy differences, and hence kink 
density at low temperatures, is also 0{Ax 2 ) at leading 
order. Corrections to the the eigenstates lead to higher 
order Ace dependences. 

More generally, given any eigenvector \tp) of the con- 
tinuum Schrodinger equation, for the specific form of SH 
of pgj) , integration by parts and use of (Eg) yields M , 



<*i«w~^<*i(£ 



(35) 



Apart from the eigenvectors, there is no temperature 
dependence in (|35|). This remarkable fact immediately 
suggests the introduction of a counterterm in the lattice 
equations which, in perturbation theory, would lead to 
the cancellation of errors of order (Acc) 2 . Modifying the 
potential as follows 



v(tf>) = v(<p) 



{Ax 



24 \8(t>J 



(36) 



leads to the cancellation of lattice dependences to order 
(Acc) 2 in a way that preserves the fluctuation-dissipation 
relation (taken at any temperature) and is thus suited for 
dynamics as well as thermodynamics. With Acc taken to 
be small enough, the leading error now becomes domi- 
nantly O(Acc) 4 . We note that unlike the situation for 
PDEs where one improves the lattice approximation for 
spatial derivatives, here a local counterterm produces the 
same effect. 

In the sp ecific case of a </> 4 potential, the coun- 
terterm (|36|) gives a new potential including the term 
—Ax 2 <j> 6 /24. The corrected potential is thus unbounded 
from below! In first order perturbation theory this is 
not a problem since the corresponding wave function is 
exponentially small in the pathological region of the com- 
pensated potential ||. However, if the full potential is to 
be used in a Langevin simulation it is clear that at suf- 
ficiently long times, the unboundedness of the potential 
implies the nonexistence of a stable thermal distribution. 
Fortunately, it is simple to estimate whether this prob- 
lem actually shows up in real simulations. The answer, 
as we show below, is that it is of absolutely no practical 
significance in the parameter range of interest. 

The resolution of this apparent difficulty brings us back 
to the validity of the expansion in Acc. So far we have 
implicitly assumed that ft ~ 1, so that Acc is the only 
small parameter and controls the order of the expansion. 
If on the other hand one wanted to work in a regime 
where Acc > ft, the whole expansion in Acc would have to 
be rederived in terms of an appropriate small parameter. 
In any case this latter regime would always constitute 
a poor approximation to the continuum: It is the Ising 
(disorder) limit of the field theory. 

A simple argument for why the counterterm works at 
small temperatures, meaning Acc <C ft, is the follow- 
ing. Consider a temperature large relative to the po- 
tential barrier between the minima. Then, from the un- 
corrected eigenvalue equation, (0(cc) 2 ) ~ /3 _1 . On the 
other hand the value of 4>{x) 2 for which a fluctuation can 
probe the effect of the negative (j> 6 term at large (f> is 
4> 2 (x) ~ 6/(Acc) 2 . Therefore the condition for the nega- 
tive <fi 6 term not to affect the evolution is A 2 cc -C 6/3. At 
lower temperatures it is more appropriate to explicitly 
calculate the Kramers escape rate [|l7j across the bar- 
rier separating the metastable and unstable regions of 



5 



the compensated potential. Assuming the lattice sites 
to be uncoupled (this gives an overestimate of the true 
rate), the calculation yields Tk ~ exp(— 4(3/3(Ax) 4 ), 
which turns out to be vanishingly small in practice: For 
Ax = .5, P = 5, and a lattice size of 10 6 points, the 
probability of an escape at a single site per unit time 
is only ~ 10 -41 . In our numerical calculations we have 
verified that the counterterm can indeed be successfully 
used in the appropriate circumstances with no hint of 
any instabilities. 



IV. THE ONE LOOP APPROACH 

In contrast to the above considerations, the 1-loop 
counterterm proposed in Ref. Q arises from the conjec- 
ture that the leading dependence of the partition func- 
tion on Ax coincides with the most divergent term for the 
same theory in higher dimensions. Although the relevant 
computations are well known we will spell out some of 
the steps to make every assumption clear. The basic idea 
is to start again with the canonical partition function: 



Z = N D<j>e 



-ps[<t>] 



(37) 



The field (j> is then decomposed into a background field 
fa and a fluctuation field Xi 4> — fa + X-, an d assuming 
the fluctuations to be small, expanded around fa: 



S[fa + X ]~S[fa} + 



ss 



' \x=o 



x - 



1 PS_ 

2 X ~5fai X 

T \x=o 



(38) 



If fa is an extremum of S[<f>] then the first term van- 
ishes. Under this assumption 



Ne -0S[4> b ] I Dxe ^ , X=0 X 



(39) 



Because 



s J s 



'lx=c 



is independent of x the functional inte- 
gration is strictly Gaussian and can be performed exactly: 

N [ Dxe~ P * x %? lx=0X = Det" 



4 I • (40) 



Here we have adopted the usual normalization to the free 
theory. The action S = Sq + Si, was decomposed into the 
action for the free theory So (gradient and mass terms) 
and the interactions Si. Primes denote functional deriva- 
tives relative to <j). 

This can be written as 



Dct 



-1/2 



S 



^ =Dct- 1 / 2 (I + K) 



= e -iTrlog(l+if) 



(41) 



where K = Sj/S . Performing the ID fc-space trace 
integral, (m = 0), with an ultraviolet cutoff A = ir/Ax, 
we obtain one loop corrections to the potential 



7 / 7 At // 

Vil (4>, A) = V a + j^JS'^fa) - -^s-Sj (fa) . (42) 

The partition function is now approximately given by 



Z = 



-ps[4> b ]+^ ,/iPfohj- ^s' T (<fb) 



(43) 



Equations ( f42| ) and ( p3[ ) constitute the basis for the 
proposal of Ref. . In order to cancel the leading Ax 
dependence arising in this scheme, the original bare po- 
tential is modified by the addition of the last term in 
Eq. (f42|), (with a positive sign). 

Notice that while a careful accounting of the dynamics 
on the lattice yields a leading correction of order (Ax) 2 , 
regardless of any assumptions about the dominant ther- 
modynamic field configurations, the one loop procedure 
leads to a correction of order Ax. In contrast to the 
correct answer discussed in the previous section, the one 
loop procedure gives no corrections for the free theory 
since in this case Si = 0. 



V. COMPARISON WITH NUMERICAL 
SOLUTIONS 

Accurate Langevin studies of even one-dimensional 
field theories require large lattices and long running 
times. It has only recently been realized, by compari- 
son against exact analytic results for nonlinear field the- 
ories, that fairly large errors (e.g., 30% or greater in the 
kink density) can easily arise if numerical studies are not 
carried out with careful error control methodologies [|ll| . 

In order to test the predictions of the previous sections, 
we ran large scale Langevin evolutions with typical lattice 
sizes = 10 6 , and with the time step related to the 
lattice spacing by At = 0.05 Ax 2 . 

A first test which allows comparison against exact an- 
alytical results are the lattice dependences for the linear 
SPDE (free theory) defined by V(<j>) = (j) 2 /2. Figure | 
shows the plot of the thermal equilibrium 1 — (</> 2 ) versus 
Ax. The numerical data are in excellent agreement with 
the (exact) theoretical predictions. 

In the more general case of a nonlinear SPDE, we can- 
not expect explicit exact solutions for arbitrary Ax, but 
thermodynamic quantities can be obtained to order Ax 2 
from the eigenvalues of the perturbed Schrodinger equa- 
tion extracted from the transfer integral, as described in 
Section III. In the case of predictions for the kink den- 
sity, precise comparison with numerical results has not 
been possible until recently, partly due to the difficulty 
of counting the number of kinks in a noisy configuration. 
The correlation length is, however, a well-defined quan- 
tity at any temperature, independent of kink-counting 
schemes. 

We extract the correlation length Aoo from the numer- 
ically determined field configurations as follows. Let 



c(iAx) = {<j>(j)(j>(j + i)), 



(44) 



G 



and 



X(x) = Ax log 



c(x) 



^c{x + Ax) 
The correlation length is lim^^oo A(x): 
{<p(Q)<f>(x)) -> exp {-x/Xoo) , 
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oo. 



(45) 



(46) 



The correlation function c(x) is in general a sum of ex- 
ponentials (the smallest exponent being the correlation 
length). For values of x much smaller than the correla- 
tion length, therefore, X(x) < A^. In practice, for large 
x, finite statistics mean that the ratio in (^5|) cannot be 
evaluated precisely. One therefore evaluates the correla- 
tion length by plotting X(x) versus x and looking for a 
plateau at intermediate values of x. 
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FIG. 1. The dependence of 1 - (4> 2 ) on Ax, for 
V(4>) — <f 2 /2. The numerical results (•) at /3 = 2 are com- 
pared with the exact equilibrium result (|3l|) (solid line) . The 
dashed line shows the Taylor expansion of (pw to order Ax 2 . 
Statistical error bars are not shown if they are smaller than 
the symbol size. 

We measured the correlation length using three dif- 
ferent Langevin evolutions: (a) A standard simulation 
using a second-order stochastic Runge-Kutta integrator; 
(b) A simulation with the counterterm (|36j); (c) A simula- 
tion with the counterterm proposed in Ref. |12|. Results 
for Ax = 0.5 are shown in Figure |^. The counterterm 
( |36| ) shifts the result from the Langevin evolutions (a) 
on the lattice very close to the exact continuum result, 
shown as a dashed line. The standard simulation over- 
estimates Aoo, whereas the one loop counterterm results 
in an underestimate with an error larger than the "bare" 
simulation (a) without any counterterm. 



1 7 - 



1 6 



1 5 




1 



1 5 

x 



20 



25 



30 



FIG. 2. The correlation length A computed with the bare 
potential (•), the counterterm of ( |36[ ) (o) and the one loop 
counterterm (A), for Aa; = 0.5 and /3 = 5. The dashed line 
shows the continuum exact result, computed via the transfer 
integral. 
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FIG. 3. Early evolution of the space-averaged mean value 
of (f> for different values of the lattice spacing Ax. From top to 
bottom, the lattice spacings are: Ax — 0.25, 0.5, 0.75, 1.0. 
We used lattices of 1048576 points and At = 0.05Ai 2 . 
/3 = 10, 77 = 1. 

As a further test we repeated, with large lattices and 
smaller time steps, a numerical experiment presented in 
Ref. JT^ |. The initial condition is chosen uniform at the 
minimum of V ((/)), 4>o = ~ lj the system is then run for 
a short time (before any kinks appear) so as to observe 
the relaxation to a mean value <p m . Although this does 
not result in a strictly thermalized configuration, small 
wave-length fluctuations quickly display a thermal spec- 
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trum. (In other words, the "phonon" relaxation time is 
much smaller than the timescale for kink nucleation.) In 
Figure [| we show the value of (0) as a function of time 
for four values of Ax. From the plateau for moderate 
times, we can obtain a fairly precise estimate of <p m . As 
a cautionary note, we point out that at small Ax, a small 
stepsize is also needed (see Figure ^) . 

It is possible to employ a Gaussian approximation (fol- 
lowing Ref. JlO| ) to obtain a rather good estimate of <f> m 
as a function of Ax, the result being shown in Figure |^: 
The leading dependence, both analytically and numeri- 
cally, is clearly quadratic in Ax. To obtain the analytic 
result, we use the fact that the probability density of 
4> is the square of the ground state of the Schrodinger 
equation (|2q). (This density emerges from dynamic sim- 
ulations or calculations; it is not an input to numerics.) 
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FIG. 4. Early evolution of the space-averaged mean value 
of 4> f° r two values of the time step At, with (5 = 10, tj = 1. 

We proceed further by using a Gaussian ansatz [jlO| for 
the ground state eigenfunction: 



,.,„(,,) = ( "Yexp( i<>.:-, G:,r;.. 



(47) 



The parameters il and 4>o are obtained by minimizing the 
energy 

/oo 
4,%(<t>)H(<t>)d<p. (48) 
-OO 

For large (3 the two free parameters are related by 

= /3(3^-l) 4 +O(l), (49) 

and 

E = -l<& + 1$ + /r x ~(3^ - 1)' + o{fr 2 ). (50) 



The dependence of <fim on Ax can now be estimated us- 
ing the Gaussian approximation of ([l7|). At finite Ax, we 
replace V(4>) in @ by V((b) + ±Ax 2 (^V(4>)) 2 . Mini- 
mizing ([50]) with respect to <po gives 



<^o(Ax) = o (O)+/3- 1 Ax 5 
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MV2 



+ o(r 2 ), (5i) 



which is plotted in Figure ^|, in excellent agreement with 
the numerical results. 
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FIG. 5. Space- and time-averaged mean value of <f> for 
different values of the lattice spacing Ax. The solid line is 
the large-/3 estimate ( pi] ) obtained from the Gaussian ansatz 
( p7| ) with P — 10, rj = 1. Statistical error bars are not shown 
if they are smaller than the symbol size. 



VI. CONCLUSIONS 

We have presented a complete procedure to identify 
space and time discreteness effects in Langevin studies 
of 1+1-dimensional field theories on the lattice. This 
scheme permits the determination of the correct con- 
tinuum limit of the theory in thermal equilibrium. In 
particular, we have shown that for the standard spatial 
discretization of the Langevin equation, quantities of in- 
terest such as the kink density and the expectation value 
of the field and its variance differ from the continuum 
values by terms of order Ax 2 . High resolution numeri- 
cal results are in excellent agreement with our analytical 
predictions. 

In any numerical Langevin evolution errors result from 
the necessary discretization of a field theory in both time 
and space. The effect of the former is to modify the form 
of the canonical distribution as seen from the solution 
of the corresponding Fokker-Planck equation. The use 
of higher order timestepping algorithms can render this 
error subdominant when compared to errors arising from 
the discretization of the spatial lattice. 
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This spatial discretization error can be computed sys- 
tematically in powers of (Ax) 2 via the use of the transfer 
integral to solve for the partition function on the lat- 
tice. This procedure leads to the identification of a sim- 
ple local counterterm which in turn permits the practical 
elimination of the leading order lattice error in Langevin 
evolutions at low temperature. 

For the <p A theory in one space dimension, the density 
of kinks converges to a well-defined value at any temper- 
ature low enough that kinks are clearly separated from 
small wave-length fluctuations (or "phonons" ) . In prac- 
tice this is essentially the range of temperatures where 
the dilute gas approximation (which is equivalent to a 
WKB solution of the transfer integral) is valid. Preci- 
sion calculations over a wide range of temperatures that 
agree with transfer integral predictions are reported in 
Rcf . [prf . For quantities that are defined unambiguously 
at arbitrary temperatures, such as the correlation length, 
results based on the WKB approximation to the transfer 
integral will fail at sufficiently high temperatures. This 
is independent of lattice errors and does not preclude an- 
alytical and numerical study using lattice simulations. 

Our results disagree with those reported by Gleiser and 
Muller based on a one loop counterterm A critical 
examination of their proposal has shown that it does not 
in fact constitute a scheme for the control and elimina- 
tion of lattice errors. We have also carried out a direct 
comparison of our numerical results with those presented 
m Ref. @ and are led to conclude that their data must 
have been of insufficient quality to quantitatively char- 
acterize lattice spacing dependences. 

Finally we wish to note that our methods strongly de- 
pend on the use of the transfer integral to solve (exactly) 
for the non-perturbative field thermodynamics. The ap- 
plication of the procedure described above is in general 
guaranteed for any (local) field theory in one spatial di- 
mension, requiring only the choice of the appropriate po- 
tential. 

In higher dimensions such a solution becomes increas- 
ingly difficult. Nevertheless, for the particular case of two 
dimensions several lattice models can be solved exactly, 
precisely by applying the transfer integral technique [jl8| . 
The eigenvalue problem, elegantly posed in terms of an 
ordinary time-independent Schodinger equation in one 
dimension, now amounts to solving for the eigenvalues 
of an infinite matrix. Regardless of this apparent dif- 
ficulty, exact solutions are known in several interesting 
cases, most notably perhaps for the 2D Ising model on 
a square lattice. It is therefore conceivable that the de- 
tailed approach to the continuum in these models can be 
understood via the procedure described above. 
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